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ABSTRACT 



Context. The gas kinetic temperature in the centres of starless, high-density cores is predicted to fall as low as 5-6 K. The temperature 
gradient, which affects the dynamics and chemistry of these objects, should be discernible with radio interferometers reaching a spatial 
resolution of 1000 AU or better 

Aims. The aim of this study was to determine the kinetic temperature distribution in the low-mass prestellar core Oph D where 
previous observations suggest a very low central temperature. 

Methods. The densest part of the Oph D core was mapped in the NH3(1,1) and (2,2) inversion lines using the Very Large Array (VLA). 
The physical quantities were derived from the observed spectra by fitting the hyperfine structure of the lines, and subsequently the 
temperature distribution of Oph D was calculated using the standard rotational temperature techniques. A physical model of the cores 
was constructed, and the simulated spectra after radiative transfer calculations with a 3D Monte Carlo code were compared with the 
observed spectra. Temperature, density, and ammonia abundance of the core model were tuned until a satisfactory match with the 
observation was obtained. 

Results. The high resolution of the interferometric data reveals that the southern part of Oph D comprises of two small cores in 
consistence with the 1.3mm dust continuum map of Motte et al. (1998, A&A, 336, 150). The gas kinetic temperatures, as derived 
from ammonia towards the centres of the southern and northern core are 7.4 and 8.9 K, respectively. These values represent line-of- 
sight averages using the LTE assumption. A model using modified Bonnor-Ebert spheres, in which the temperature decreases to 6. 1 
K and 8.9 K in the centres of southern and northern core, matched the observed values satisfactorily. The southern core, which has 
more steep temperature gradient, has central density of =4x10* cm"', and the data suggests depletion of ammonia within 700 AU 
from the centre. The northern core, which is almost isothermal, seems to be less dense. The radial velocity gradients in these cores are 
almost opposite in direction, which may be an indication that turbulent fragmentation has a role in the formation of these cores. The 
observed masses of the cores are only ~ 0.2 Mq. Their potential collapse could lead to formation of brown dwarfs or low-mass stars. 

Key words. Stars: formation - ISM: clouds - ISM: individual objects: Oph D ~ ISM: molecules - Techniques: interferometric 



1. Introduction 

Dense cores within molecular clouds are birth places of stars as 
some of them are gravitationally bound and will subsequently 
collapse and form a star. Observations of starless, prestellar 
cores promote our understanding of the initial stages of star for- 
mation. Molecular line observations provide information about 
e.g. gas density, cloud structure, and kinematics of the gas. 
Furthermore, they can be utilised to probe gas kinetic temper- 
ature. Radiative transfer calculations suggest, that the external 
heating by the interstellar radiation field is attenuated by dust in 
dense cores thus causing the gas temperature to fall tow ards the 
centre of these cores (see , e.g.|Zucconi et al. 2001; Evan s et aP 
l2001l;lKeto & Field '2005'). The temperature gradient is most pro- 
nounced within a radius of a few thousand AU from the core 
centre. Far infrared- and submillimeter dust continuum studies 
have shown dust temperature to decrease towards the centre of 
cores (see, e.g. Ward-Thompson et al. 2002; Juvela et al. 2002; 
IPagani et al.l2003ll2004tlSchnee & Goodman 2005). High angu- 
lar resolution observations showing this behaviour for gas have 
not b een extensively conducted (one of the few is ICrapsi et aP 
I2007h . 



* This publication is based on data acquired with the Very Large 
Array operated by The National Radio Astronomy Observatory 
(NRAO). NRAO is a facility of the National Science Foundation opera- 
ted under cooperative agreement by Associated Universities, Inc. 



An isolated, roundish starless core with high central density 
(«c > 10^ cm-^) is a preferred choice to look for the predicted 
temperature gradient. One nearby core that fulfills these criteria 
is Ophiuchus D (LI 696 A), north from the p-Ophiuchus molec- 
ular cloud. We adopt the new dist ance estimate of 120 p c to the 
p-Ophiuchus cloud for this core (iLombardi et al.ll20p8[). There 
is no embedded protostar in southern part of Oph D ( Kirk et all 
2007). Our reference position RA = 16''28™28';9, Dec = -24°19 
'09"(J2000), coincides with the dus t emiss ion maximum of the 
850 yum SCUBA map of iKirk et a l.'' ('2005' ), the dust absorption 
maximum of 7 fim ISOCAM map of Bacm ami et al.! (2000). and 
the N2D^ line emi ssion peak according to the observation by 
ICrapsi et alJ (l2005l) . Direct evidence of very cold gas in molec- 
ular clouds was obtained bv lHariu et alJ ([2008) when observing 
the ortho-li2D^ line towards Oph D with APEX, using the same 
antenna pointing position as our reference point. They reported 
the kinetic temperature of the core derived from HaD^ linewidth 
to be 6.0+1.4 K. 

Nitrogen bearing molecular species, such as NH3, are found 
in high density cores (hc > 10^ cm"^) where most carbon bearing 
speci e s have already depleted from gas phase (e .g., Willacv et al 
'1998'; Aikawa et alj 1200 it iTafalla et al.l l2002t |j0rgensen et al 
2004). In addition to tra cing the denses t gas NH3 is also a good 
thermometer of the gas dHo et alj|1979 ^. and therefore appropri- 
ate choice for probing the temperature gradient. 
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2. Observations and data reduction 

Maps of the /:iflra-NH3 iJ,K) = (1,1) and (2,2) transitions at fre- 
quencies 23 694.4955 MHz and 23 722.6336 MHz were made 
with the Very Large Array (VLA) in five observing slots, each 
six hours in duration, between 9th and 15th June 2008 in the ar- 
ray's DnC configuration. The synthesised beam of VLA in this 
hybrid compact configuration was 3'.'72 x 2'.'35 (position angle 
PA = 37°.2 east of north) with natural weighting applied. The 
FWHM of the primary beam is 2' at this frequency. In order to 
improve the signal-to-noise ratio of the (2,2) map we used a cir- 
cular clean beam with FWHM of 10" to construct the clean im- 
age. For consistency we also cleaned the (1,1) map with the same 
beam to ensure identical resolution of the maps for the analysis. 
The transitions were observed simultaneously in 21F mode with 
128 channels and a bandwidth of 1.5625 MHz in both spectral 
windows (resulting in a channel width of 12.2 kHz or 0.154 km 
s~'). The bandwidth was enough to include the main hyperfine 
group and the inner Fi - 1 - 2 and Fi -2-1 satellite groups 
of the (1,1) line. 

Temporal variations in gain and phase were calibrated 
every 15 minutes using quasar 1626-298 (Fi.3cm = 2.57 J>[3). 
Bandpass- and absolute flux calibration were carried out in 
the beginning and the end of each 6 hour period using quasar 
1331H-305 {Fi,2cm = 2.59 Jy). The uncertainty of flux calibration 
was estimated to be of the order of 15%. At the time of the obser- 
vation some of the antennas of the array were replaced by new 
Expanded VLA antennas, for which reason the online Doppler 
tracking was not available. In order to avoid abrupt phase jumps 
the sky frequency was calculated by hand for each phase cali- 
brator - source - phase calibrator triplet. 

The visibilities were calibrated and data reduction and 
imaging was carried out with Common Astronomy Software 
AppHcations (CASAfl The noise rms of the data cube cleaned 
with the 10" beam is 1 mJy/beam, which corresponds to about 
22 mK. We were unable to estimate the missing flux since no 
single dish data was available (the maximum in the full resolu- 
tion clean map is 0.11 K km s ' whereas the negative contours 
reach a level of -0.06 K km s '). This should not affect the uncer- 
tainties of physical quantities derived from line intensity ratios, 
e.g. optical thickness and kinetic temperature, unless there is dif- 
ferential flux loss between (1,1) and (2,2) spectral lines. 

We fitted the hyperfine structure of NH3(1,1) -line using our 

-minimisation routine to obtain estimates of line total optical 
thickness ri i , velocity i^lsr, and linewidth (FWHM) Av. The hy- 
perfine structure of (2,2) -line could not be detected, so the main 
component was fitted with a single Gaussian line profile. 

3. Results 

3. 1 . Integrated intensity maps 

The integrated intensity maps are shown in Figs.[T]and|2] The full 
angular resolution map of the ammonia (1,1) transition is shown 
in Fig. [T] The signal-to-noise ratio of the (2,2) transition was 
very small when cleaned to the full angular resolution, and there- 
fore its map is not shown. In Fig.|2]the integrated intensity maps 
of (1,1) and (2,2) transitions smoothed to 10" are presented. The 
emission of (1,1) line was integrated over all available hyperfine 



VLA calibration manual: 
|http : //www . via . nrao . edu/astro/calib/manual/csource ■ htinl| 

^ CASA is being developed in a collaboration led by NRAO as 
a main post-processing package for ALMA and Expanded VLA, see 
[httpV/casa. nrao.edu, 




60 40 20 -20 -40 
RA offset (") 

Fig. 1. Integrated intensity map of Oph D in (1,1) transition of 
NH3. The VLA synthesised beam, 3'.'72 x 2'.'35, is indicated in 
the bottom right corner The contours are from 20% to 90% with 
steps of 10% of the map peak value 0. 1 1 K km s"' . The reference 
position is RA = 16''28'"28!9, Dec = -24° 19 '09"(J2000). 

Table 1. Velocity gradients of Oph D measured from the 
NH3(1,1) lines. 



Source" 








[km s"' pc"'] 


n 


N 


6.71+0.01 


178.7+0.4 


S 


2.10+0.06 


55.8±1.5 



Notes. The northern peak is indicated with N, and southern with S. 
Direction of increasing velocity measured west of north. 

components (main group and the inner satellite groups) and (2,2) 
line over the main component. The RA and DEC offsets in these 
maps are relative to our reference position. 

The integrated intensity maps show two prominent lo- 
cal emission maxima within the mapped area, especially in 
NH3(1,1) maps. Hereafter, we will refer to these as the north- 
ern and the southern peak. Such structure has not been seen in 
previous submillimeter or infrared observations no r in molecu- 
lar line mappings of this core, with the exception of iMotte et al.l 
119981) who report two cores in the southern end of Oph D in 
their 1.3 mm continuum map (D-MMl and D-MM2, see their 
Fig. 2). Typically the southern end of Oph D has been characte- 
rised as a single prestellar core. The structure will be discussed 
in more detail in section 4. 

3.2. Kinematics of tiie gas 

The kinematics of the gas was studied with the aid of Gaussian 
fits to the hyperfine structure of the NH3(1,1) line which pro- 
vide accurate estimates of radial velocity and line widths. The 
resulting Icr uncertainties of these parameters are typically bet- 
ter than 0.01 km s '. The hyperfine fits were used to form an- 
other image cube with a single Gaussian line profile normalised 
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Fig. 2. Integrated intensity maps of Oph D in (1,1) (left) and (2,2) (right) transitions of NH3 smoothed to 10" to improve the signal- 
to-noise ratio of (2,2) line. The 10" Gaussian clean beam is indicated in the bottom right corner The contours are from 10% to 90% 
with steps of 10% of the map peak values 3.5 K km s"' and 0.1 1 K km s"' for (1,1) and (2,2) lines, respectively. 
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Fig. 3. Channel map of Oph D in the (1,1) transition of NH3. 
Channel width is 0.154 km s 



with the integrated intensity of the full NH3(1,1) line. This im- 
age cube describes the radial velocity distribution of the ammo- 
nia. The velocity channel map based on this image cube is shown 
in Fig. [3] The channel width corresponds to the original resolu- 
tion of the VLA correlator. The map covers velocities from 3.29 
to 3.76 km s ' in four channels. The emission of the northern 
peak moves from northeast towards southwest as the velocity 
increases, whereas the southern peak emission appears to move 
slowly westwards along a roughly east-west path. 



We created a local velocity gra dient map for both tran sitions 
following the method presented in lGoodman et al.l (Il993l) . Each 
velocity gradient is determined using a nine pixel grid centred at 
the pixel under consideration. Furthermore, we required that at 
least seven of the surrounding pixels have a usable velocity esti- 
mate available for the least squares fit. The results are presented 
in Fig 3] where the velocity gradient vectors of (1,1) transition 
are superimposed on the map of integrated intensity. 

In the velocity gradient maps one can see two distinctly dif- 
ferent velocity fields on top of the southern and the northern 
peaks, which seems to suggest dissimilar rotation of gas near 
these peaks. The effect can also be seen in the map of the (2,2) 
transition, although not as prominently. The velocity gradients in 
the vicinity of both peaks seem homogeneous enough to apply 
the solid-body rotation fit also in a larger scale. The overall total 
gradient is determined for the two locally uniform areas sepa- 
rately instead of giving one value to the whole clump. The mag- 
nitude and direction of total gradients for northern and southern 
peaks are given in Table [1] 

The magnitudes of the total velocity gradients, Q - 6.71 
and @ - 2.1 km s ' pc of the northern and the southern 
peak are comparable to the values typi cally detected in oth er 
cores using nitrogen co mpound tracers (iCaselli et al.l l2002allbt 
ICrapsi et"d]|2004 l2007i) . bu t the agreement with previous ob- 
serva tions of Oph D is poor (iGoodman et al.l[T993t ICrapsi et alJ 
'2005). This is probably due to significant difference in the angu- 
lar resolution of the data, being ~ 80" (Benson & Myers 198^ 
see their Fig. 21) for NH3 data used in Good man et al.l (11993h 
and 26" for the N2H+(l-0) data used in Cra psi et al.1 (l2005h 
With these beamwidths the velocity field is averaged consider- 
ably, and more importantly the reported magnitude of the total 
gradient of the whole core will not be a good figure of rotation 
for this object, since, as we have shown, there exists two distinct 
velocity regions with almost opposite gradient directions. 
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Table 2. NH3(1,1) and (2,2) fits and derived parameters in the northern and the southern local column density maxima of Oph D 
using the 10" resolution map. 



Source / Line" ti^sR A« fZ 7km MNHs) 

[K\ [kms-'] [kms-'] [K\ [K] [K] [10'^ cm'^] 

N,NH3(1,1) 3.52+0.06 3.39±0.001 0.18+0.003 10.37±0.39 6.05+0.50 8.59+0.37 8.87+0.56 9.60+2.97 

N, NH3(2,2) 0.28+0.10 3.38+0.02 0.27+0.05 

S, NH3(1,1) 3.33+0.12 3.47+0.002 0.16+0.005 7.16+0.52 5.70+0.50 7.32+0.59 7.43+1.12 9.44+2.75 

S, NH3(2,2) 0.12+0.10 3.48+0.03 0.25+0.09 



Notes. The northern peak is indicated with N, and southern with S. The exact positions of the peaks are: RA = 16''28"29:3, Dec = -24° 19 
'03"(J2000) for the northern and RA = 16''28"28';8, Dec = -24° 19 '23"for the southern peak. <*> The sum of peak optical thicknesses of all (1,1) 
hyperfine components. This is twice the 'main group opacity' used in many papers. 
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Fig. 4. Velocity gradient map of NH3(1,1) transition of Oph D. 
The arrows depict the direction (of increasing radial velocity and 
the magnitude of the velocity gradients. An arrow with a length 
of 10" represents a gradient of 14.3 km s"^ pc 



3.3. Ammonia column densities and tine kinetic temperatures 
using tiie LTE assumption 

We fitted the hyperfine structure of NH3(1,1) line using ax^ rou- 
tine, and the (2,2) line with a simple Gaussian profile. Fits were 
performed for pixels having signal larger than 3cr„oise level. We 
made the usual assumption, that the different hyperfine states 
are populated according to LTE. We derived the excitation tem- 
perature Tex, rotational temperature T12 and the total (ortho -H 
para) ammonia column density A^(NH3) following the proce- 
dure described by Hariu et al. ( 1993). In the calculation of to- 
tal ammonia column density we assumed the NH3 ortho-to-para 
ratio to follow the LTE conditions. The kinetic temperature of 
the gas was calculated from the rotational temperat ure using the 
semi-e mpirically derived analytic expression from lTafalla et al.l 
(12004 



l-(ri2/42)ln[l + l.lexp(-16/ri2)]' 

which is based on Mont e Carlo simulations using NH3 - H2 col- 
lision coefficients from iDanbv et all (119881) and is found to be 




Fig. 5. Spectra (gray bars) and chi-square fits (black solid line) 
of NH3(1,1) and (2,2) from the centre of the northern core (from 
the 10" resolution map). The (2,2) spectrum is off'set by -1 K for 
clarity. The splitting of the (1,1) components is due to hyperfine 
structure, and is not to be taken as a sign of self-absorption. 



suitable for d ense cores h aving temperatures less than 20 K. 
According to iMaret et aP (12009) . who present new coUisional 
rates for NH3 - H2, their new and Danby's coefficients agree 
within a factor of two. They also found, that non-LTE Monte 
Carlo radiative transfer simulations using their new and Danby's 
coeflicients produce practically equal results for a parameter set 
typical of a dense and cold prestellar core. Thus, Eq. ([T]i charac- 
terises the relation between rotational and kinetic temperature of 
the gas robustly. 

The results of the hyperfine fits and the physical parameters 
derived for the northern peak (RA = 16*'28™29';3, Dec = -24°19 
'03") and the southern peak (RA = 16'>28™28^;8, Dec = -24°19 
'23") peaks are presented in Tabled The maps of NH3 optical 
thickness t and column density A^(NH3) are shown in Fig.|6l 

3.4. Mass and fractionai abundance 

In the abse nce of continuum data of our own, we use the 1.3 
mm data o f lMotte etal.l ( [T998h to estimate masses and fractional 
abundance of ammonia of our sources. Despite the difference of 
5" between the coordinat es of our core cent res and the D-MMl 
and D-MM2 reported by iMotte et all (1 19981) . we consider them 
the same cores. 
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If we assume the dust is isothermal (having temperature To) 
and the background emission has been subtracted from the ob- 
served intensity, the intensity of dust emission /J?"'*' can be calcu- 
lated using the equation of radiative transfer: 



(2) 



Thermal dust emission in the mm/submm region is optically 
thin, so Eq. (|2]i becomes /jj"'*' x By{Ti)Ty. The optical thickness 
Tv is given by the product of total surface density 2, the mass 
absorption coefficient per unit mass of dust /cj and the dust-to- 
gas mass ratio R^. The surface density can be expressed with the 
aid of column density S = A^(H2)/iH2'WH, where //h, is the mean 
molecular weight per H2 molecule and oth is the mass of the hy- 
drogen atom. The column density as a function of dust intensity 
and temperature becomes 



rdust 



A^(H2) 



By{Ti)^iYi2fn\^KiRi 



(3) 



The 1. 3 mm flux density in a 15" beam reported bv lMotte et alj 
(119981) is 45 +5 and 45 ± 1 5 mJy beam ^ for southern and northern 
core, respectively. We assume that the gas and dust temperatures 
are equal, and use an average value of Td = 8 K for the dust 
tem perature for both cores. Ado pting values of - 0.5 cm^ 
g~' dOssenkopf & Hennin3ll994^ and - 0.01, the molecular 
hydrogen column density for the cores becomes 



A^(H2) = 5.3 : 



: 2.0 X 10^^ cm ^ for the northern and 



N(H2) 



5.3 + 0.6 X 10^^ cm ^ for the southern core. 



The fractional abundance of ammonia can be estimated using 
the calculated value of molecular hydrogen column density and 
our ammonia column density (Table|2]i, 



X(NH3) = 



A^(NH3) 
A^(H2) 



1.8 X 10" 



The uncertainty of factor ~2 is typically associated with the 
dust opacity coefficient. This affects directly the column density, 
ma ss, and abun d ancg e stimates. 

iMotte et all (Ii998b give mass estimates for both cores based 
on integrated flux densities. The reported masses of the cores 



are 0.15 Mq and 0.16 Mq for the southern and northern cores. 
Based on peak flux densities, the mass within the 15" beam can 
be estimated with 

M=^^, (4) 

where d - 120 pc is the distance to the source. This gives a mass 
estimate of M = 0.10 M0 for both cores. The 15" beam, how- 
ever, encloses the cores only partially, and therefore the mass 
estimate is probably too small to describe the whole cores. 



4. Modelling 

4.1. Description of the model: BE splieres 

The southern tip of Oph D cloud has in general been considered 
a single prestellar core. Our VLA data, however, indicates the 
existence of two cores instead of one. This hypothesis is backed 
up by the velocity gradient map (Fig. |4| as well as by the two 
separate local maxima of optical thickness and column density 
seen in Fig. |6] The velocity gradients seem to imply, that the 
direction and magnitude of rotation of the gas within the area 
of these two hypothetical cores is markedly different from each 
other The secondary peaks of optical thickness and column den- 
sity south from the m ost dominant peak also point towards the 
presence of two cores. Steinacke r et alj ([2005) modelled Oph D 
cloud with a complex 3D substructure based on ISOCAM 7 and 
15 jum maps and a 1.3 mm map from IRAM 30 m telescope. 
In their model the southern part of the cloud also contains two 
dominant density concentrations. Also Motte et al. ( 1998) report 
two cores (D-MMl and D-MM2 in their Fig. 2 and Table 2) in 
the southern part of Oph D. 

We attempt to model the southern part of Oph D with two 
cores, each in hydrostatic equilibrium and with a radial tem- 

Ierature gradient, the so called modified Bonnor-Ebert spheres 
Evans et all I2OOII: IZucconi et all l200ib with radii 2200 AU 
(northern core) and 2100 AU (southern core). The central posi- 
tions of the core models were chosen to coincide with the lo- 
cal kinetic temperature minima derived from the observation. 
The calculation of the gas density and temperature profiles for 
both core models was done iteratively until satisfactory conver- 
gence was achieved. The initial core model was calculated by 
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Fig. 7. Line-of-sight average profiles of column density, kinetic 
temperature and optical thickness (rtot) towards three radial di- 
rections (E of N) from northern core centre outwards. The values 
derived from the observation are marked with small gray circles 
and from the model with a solid line. 



Fig. 8. Line-of-sight average profiles of column density, kinetic 
temperature and optical thickness (rtot) towards three radial di- 
rections (E of N) from southern core centre outwards. The values 
derived from the observation are marked with small gray circles 
and from the model with a solid line. 



ISipila et al ] (1201(1 see Sec. 2.1. therein), who starting from an 
isothermal sphere derived the dust density and temperature pro- 
file for a core with a radius of 2400 AU. The southern core den- 
sity is high enough > 10^ cm'-^) to assume the d ust and gas 
temperatures to be equal (iBurke & Hollenbachll983h . but for the 
northern core («c ~ 10^ cm"^) the situation may not be as clear. 
Despite this uncertainty we used Sipila's dust temperature and 
density profile, scaled to the radii of our core models, as a first 
estimate of our cores. 

In order to derive the gas density profile the f ollowing 
mod ification to the original substitutions introduced bv iBonnon 
11956) was made to include the non-isotherm al temperature dis- 
tribut ion T(r) of the core in the calculations dSipila et al.ll2010l 
1201 Ih : 

/ T \ 

(5) 



(6) 



p(r) = A I e 

^ \Tir)! 



AnGmA 



where Tc and A are the temperature and density at the centre. The 
resulting equation of hydrostatic equilibrium was solved numeri - 
cally using boundary conditions presented in'Sipila et al.l (l2010l) . 
and the obtained density and temperature profiles were fed to the 
radiative transfer Monte Carlo code (Juvela|[1993) to simulate the 
spectra. We use molecula r data provided by t he Leiden Atomic 
and Mole cular Database (Schoier et al.ll2005h . The method de- 
scribed in'Ke to et al. liooi is used for handling overlapping hy- 
perfine emission lines. For the simulation we constructed a 3D 
model which consists of these two cores with a projected separa- 
tion of 2600 AU in the plane of the sky. In addition to the density 
and temperature, the fractional abundance of ammonia X(NH3) 
in our model can be adjusted, as well as the turbulent velocity 
within the core. The non-thermal velocity dispersion component 
AiiNT is calculated from the observed linewidths A^obs with an 
equation: 



A^NT = ^ (Auobs)^ -81n2 



(7) 



The prevailing non-thermal velocity dispersion in the northern 
and southern cores is between 40-50 m s"' (the isothermal sound 
speed for a core with average gas kinetic temperature of 8 K is 
Cs ~ 170 m s"'). The turbulent velocities in our simulation model 
are therefore also kept at a low value. 

4.1.1. Comparison between the model and observation 

In dense prestellar cores the fractional abun dance of NH^ may 
increase towards the core centre (see e.g. iHotzel etal.M2004 
and references therein). In more evolved cores, however, when 
the central density increases above ~ 10^_cnr_^, also am- 
monia starts to vanish from the gas pha se dAikawa et al. 20051 
iFlower et al.l2006l;lMiettinen et al.l2010h . A variety of fractional 
abundance values was tested in core models, using the estimate 
X(NH3) w 1.5 X 10"^ discussed in section 3.4 as a starting point. 
In each simulation the abundance radial profile shape followed 
the current density profile shape, except in the case of depletion 
which approximated as a steep descending slope. 

The best-fit parameters of the northern core turned it into 
a nearly isothermal sphere; the temperature at the centre is 8.8 
K and increases only 0.4 K towards the outer edge. The op- 
tical thickness and column density could not be fitted simul- 
taneously to match the observed values at the centre of the 
cores, and a relative error of ~ 15% remains. The central den- 
sity of the model is «c =2x10^ cm^^ and the fractional abun- 
dance gradually increases towards the centre peaking at the value 
X(NH3) = 2 X 10"^. Mass of the northern core model is only 
0.045 Mq, and could not be increased while maintaining a rea- 
sonable agreement with the observation. 

In order to model the low temperatures observed in the 
southern core the model central densit y had to be i ncrea sed to 
4 X 10'' cm \ similar to the results of Hariu etlD (l2008i) . The 
best-fit temperature profile has a temperature decrease from 10 
K to 6.1 K from the surface to the centre. The fractional abun- 
dance of ammonia was then used to control the optical thickness 
and column density. It was quickly noticed that a swift decrease 
in abundance within the innermost 700 AU was required to ob- 
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Distance [AU] Distance [AU] 

Fig. 9. The density (solid line) and temperature (dash line) pro- 
files of the southern (left) and northern core (right) BE model. 



tain a match with the observed values. The abundance increases 
to a value of ^(NHa) = 2.5 x 10"' before the depletion starts to 
take effect. The mass of the southern core model is 0. 16 Mq. 

The structure of the southern part of Oph D is certainly more 
complex than what can be modelled using two spherical cores in 
hydrostatic equilibrium. Nevertheless, a relatively close agree- 
ment with the observation was achieved with this simple model. 
The spectra produced by the radiative transfer calculations were 
analysed in the same way as the observed spectra. The results are 
presented in Figs. Q and [8] In both figures the column density, 
kinetic temperature, and optical thickness are plotted in three 
radial directions from the core centre towards the edge. Each 
subplot contains the observed values in a sector +20° from the 
mean direction plotted with small gray circles together with their 
Icr uncertainties. The corresponding simulation result is shown 
with a solid line. The temperature and density of northern and 
southern core models are presented in Fig.|9l 

The 3 D model was also u sed to reproduce the H2D^ line ob- 
served bv lHariu etaP (l2008l) w ith APEX. The fra ctional abun- 
dance of H2D^ was taken from ISipila et aTl ( 1201 Oh and applied 
for both cores. Otherwise the model was exactly the same as 
in NH3 simulations. A single pointing spectrum at our and 
Harju's reference position using APEX antenna beam FWHM 
of 17" was simulated, and a linewidth of Au - 0.35 km s"' 
was obtained. This is exactly the sa me l inewidth produced 
by a single hydro static core model of fHariu etal.1 fcOOSh and 
ISipila et aU (I2OIOI) in their respective simulations of this obser- 
vation. Although the observed narrow linewidth Au = 0.26 km 
s~' could not be reproduced with our two core model, we can at 
least see that the simulation is in line with the previous single 
core models, and did not introduce further deviation from the 
observed value. 

4.2. Description of an alternative model: spheroid cores 

The modified Bonnor-Ebert sphere models reproduce the ob- 
served spectra reasonably well, and is also an apparent choice 
for the roundish southern core. The shape of the northern core, 
however, appears more elongated in the VLA map. An oblate 
spheroid, whose axis of rotation is on a plane perpendicular to 
our line of sight, would describe the shape and velocity field 
of the northern core perhaps better than a sphere. To test what 
effect the rotation of the cloud has on its shape, we run simu- 



lations with the Pencil using a 3D-grid. Pencil Code is a 

modular solver of magnetohydrodynamic equations using a high 
order finite difference scheme. The solution method enables a 
straightforward combination of various "physics modules", like 
self-gravity, and various equations of state to the general equa- 
tions describing the time-evolution of a hydrodynamic flow. As 
an initial condition, we used an isothermal Bonnor-Ebert sphere, 
where an initial rotation rate corresponded to the angular ve- 
locity of the northern core Q. = 2.2 x 10 rad s No extra 
rotational energy was put into the system during the time devel- 
opment of the simulation, and magnetic fields were neglected. 
The rotation rate foO - 0.27, where the dynamical timescale for 
BE sphere is fo - (47rGpc) '^^ (Banerjee et al. 2004), is quite 
fast, but still not enough to play alone a significant role in sup- 
porting the core against gravitational collapse. The dynamics of 
the model over time were observed, and with this angular ve- 
locity, according to calculations, the core seems to settle in an 
spheroidical equilibrium state. 

The experiment indicated that the Bonnor-Ebert sphere 
transforms into a spheroid (with an axial ratio of the projected 
ellipse less than 2: 1). The convergence in the simulation is slow, 
however, and still after 1 Myr there are oscillations. At this point 
the central density has decreased to one fifth of the original 
value, and the density profile is somewhat constant along the 
longer axis until a rather sudden drop occurs at the edge of the 
core. In the direction of the axis of rotation the shape of the ori- 
ginal hydrostatic equilibrium density distribution remains about 
the same, just rescaled. 

In this alternative model a slightly spheroidical shape was 
adopted for the southern core as well. The density and tempera- 
ture profiles of the spheroids are, as was indicated by the simula- 
tion conducted with the rotating sphere, radially scaled versions 
of those of modified Bonnor-Ebert spheres. In this model the 
central densities were = 1.0 x 10^ cm"-' and - 1.4 x 10^ 
cm"^ for the northern and southern cores. The fractional abun- 
dance in both cores has a decrease in the centre, i.e. ammonia is 
depleted. The northern core abundance reaches a peak value of 
X(NH3) = 1.1 X 10"^ 1000 AU from the centre and the southern 
core abundance ^(NHs) = 7 x lO"** 700 AU from the centre. 

The best-fit temperature profile of the northern core is 
steeper than in the previously discussed Bonnor-Ebert sphere 
model. The model gas temperature decreases to 7.4 K in the cen- 
tre of the core, and rises to 10.2 K at the surface. The southern 
core temperature profile has a minimum at the centre of 6.4 K 
and rises to 8.0 K towards the surface. Masses of the models are 
0.11 Mq and 0.07 Mq for northern and southern core, respec- 
tively. 

4.3. Comparison between the two models 

The simulation results of the two spheroid model show that the 
temperature profile of the northern core is not in good agreement 
with the observed values. The column density and optical thick- 
ness, on the other hand, are slightly better matched than with the 
modified Bonnor-Ebert sphere model. The central depletion of 
ammonia and higher gas density seem to provide more promis- 
ing match than the less dense core with incresing abundance to- 
wards the centre. 

Both models of the southern core have problems with optical 
thickness and temperature profiles. The modified Bonnor-Ebert 
sphere model seems to do better with the temperature whereas 



^ ,http : //code ■ google ■ com/p/pencil-code/| 
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Table 3. Comparison of^^ -values of fitted models. 









Til 


Northern, model 1 " 


1.38 


1.12 


1.87 


Northern, model 2 * 


1.28 


1.70 


1.61 


Southern, model 1 " 


0.91 


1.34 


1.58 


Southern, model 2 * 


1.44 


1.58 


2.03 



Notes. The model with modified Bonnor-Ebert spheres. The 
model with spheroids. 



with the optical thickness profile both models fail to explain the 

observed values consistently. 

We estimated the goodness of fit by calculating the reduced 
value between the model profile and the data: 



1 " 

-y 



(8) 



where D is the number of degrees of freedom, O the observed 
value, E the modelled value and cr^ the variance of observation. 
The -values are presented in Table |3] The comparison bet- 
ween the two models reveals no clear winner; both models have 
strong and weak properties. For the southern core, however, the 
use of a modified Bonnor-Ebert sphere leads to better results. 

It should be noted, however, that models with only one core, 
be it a spheroid or modified Bonnor-Ebert sphere, did not pro- 
duce acceptable match with the observed values. Also, if we dis- 
miss the temperature gradient of the southern core and make it 
isothermal, the temperature estimated from the simulated am- 
monia spectra is severely different from the observation. This 
indicates that a temperature gradient is necessary in the physical 
model of the southern core. The northern core is almost isother- 
mal, and its best-fit radial temperature profile has only 0.4 K 
difference between centre and surface. 



5. Discussion 

As was seen in the previous section, our distinctly different mo- 
dels both produce a reasonable match with the observation, but 
none of the two can be clearly pointed superior. The model with 
two hydrostatic cores, modified Bonnor-Ebert spheres, is, how- 
ever, on more solid ground theoretically and we will commit our- 
selves to that model in the following discussion. 

A Bonnor-Ebert sphere satisfies the following virial equilib- 
rium equation 

2r + f/ = syPext, (9) 

where Pext is the external pressure, and T and U the kinetic and 
gravitational potential energies. We neglect the effects of possi- 
ble magnetic field, and calculate the kinetic energy 

T = £therm + ^rot + ^turb = -M + -/O^ + -M- 



2 finiii 2 2 8 In 2 

(10) 

and gravitational potential energy 



U ^E, 



-GM^ 



R 



-m) 



^0 



dip 



(11) 



of the cores to examine their relations. From the relation (|9|l we 
can estimate the external pressure assuming that the cores are 
Bonnor-Ebert spheres. The centrally peaked density profile was 



taken into account individually for both cores in the numerical 
calculation of moment of inertia, /, in the equation of E^ot- The 
function /(^o) takes the density profile into account in gravi- 
tational potential energy for an isothermal core. We will use it 
although our cores have a temperature gradient (the slight dif- 
ference between the isothermal a nd non-isotherm al BE spheres 
is illustrated in Appendix B of Sioila et aD(l201 ih V The value of 
/(^o) for a homogeneous density core is 0.6. The southern one of 
our core models is slightly supercritical - 6.75 (/(^o) = 0.74) 
and the northern one has a markedly lower value = 1-49 
(/(fo) = 0.61). The energy ratios are presented in Table|4] 

The ratio of thermal energy to gravitational energy of the 
northern core is 3.89. This would suggest, that the external pres- 
sure inhibits the core from dispersing. On the other hand, when 
the core mass is estimated using data partly from the observation 
and partly from the model the energy ratio is ~ 0.5. The core 
mass in this case is estimated using the column densities within 
the area of projected ellipse in the interferometric map. When 
the average fractional abundance is estimated from the model 
Bonnor-Ebert sphere (^(NHs) = 1.8x10*^ for the northern and 
X(NH3) = 5 X 10"'' for the southern core), the total mass of the 
core can be calculated. The difference between this mass and the 
model mass for the northern core suggests that the spheroidal 
ellipsoid into which an initially rotating BE sphere evolves to, 
might be a more suitable core model. However, the mass of the 
spheroid model was also very low, only 0. 1 1 Mq. 

The modified BE sphere mass of the southern core is almost 
equal to the mass derived from the observations as described 
above. The southern core seems to be gravitationally bound, 
which can be seen from the ratio of kinetic and gravitational 
energy in Table |4] The mass of the c ore is very close to the es- 
timation given by Mot te et al.l (1 19981) . 0.15 Mq. Also the peak 
flux density of Motte's 1.3 mm continuum data was equal for 
both cores. Due to their vicinity the external pressure for both 
cores should be roughly the same. The values are close when 
the northern core mass is the one derived from the observation. 
Using the model mass the difference is significant (see Table|4]l. 
These facts further support the conclusion, that the northern core 
model is not correct, despite the fact that the match with the ob- 
served quantities is satisfactory. 

6. Summary and conclusions 

We have presented an interferometric observation of the south- 
ern part of Ophiuchus D prestellar core. We mapped the core 
with VLA using the NH3(1,1) and (2,2) inversion lines. The syn- 
thetic beam of VLA allowed spatial resolution in the order of 
400 AU, but we were forced to decrease it to about 1000 AU 
in order to improve the (2,2) line signal-to-noise ratio. From the 
observed spectra we derived the gas kinetic temperature distri- 
bution within the mapped area. The total column density of am- 
monia, the optical thickness of both lines, and the radial velocity 
distributions were also obtained. 

Our observation shows that the southern part of Oph D com- 
prises of two potentially prestellar cores. In the southern core the 
line-of-sight average gas temperature derived from the observed 
ammonia spectra drops from the outer edge towards the centre, 
reaching a minimum of Tkin = 7.4 K. Towards the centre of 
the northern core the line-of-sight average gas temperature from 
ammonia is Tkin = 8.9 K. The northern core is almost isother- 
mal with only 0.4 K increase in temperature from centre to the 
edge. The observed masses of the cores, 0.21 Mq for the north- 
ern and 0.17 Mq for the southern core, are quite small. Their 
potential collapse could lead to the formation of brown dwarfs 
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Table 4. Physical parameters of the two cores in the southern part of Oph D. 



Core 


^'-'model 




My,, 






^therm/ l^grav 1 




-£turb/|£'gravl 






[Mq] 


[Mq] 


[Mq] 


[K] 


[10* cm-^^] 








[K cm-3] 


Nortliern 


0.045 


0.21 


0.29 


8.8 


0.2 


3.89 (0.497) 


0.077 


0.045 


2.0x10" 


Southern 


0.16 


0.17 


0.21 


6.1 


4.0 


0.75 (0.413) 


1.23x10-^^ 


1.0x10-2 


9.1x10" 



Notes. The mass of the modified Bonnor-Ebert sphere used to model the core. *** The mass of the core derived from the observation with the 
aid of the model. The projection area of the core and the column density of ammonia within that area were taken from the observation. In order 
to calculate the mass, a value for the fractional abundance of ammonia was taken from the model (modified BE sphere). The average values of 
abundance over the projected area of the model were estimated to be ^(NHs) = 1.8 x 10"** for the northern and ^(NHs) = 5 x 10"' for the southern 
core. The kinetic temperature at the centre of the core model. The number density at the centre of the core model. Ratio of thermal energy 
to gravitational potential energy using the core model mass M^odeh and in parenthesis using Mobs. The external pressure, when using the core 
mass derived from observation (Mobs = 0.21 Mq), is Pcxt/^B = 9.9 x 10* K cm"^. 



or low mas s stars . The two cores were previously detected by 
iMotteet all (119981) . 

Our observation shows that there exists two rather homoge- 
neous velocity gradient fields within the mapped area. They co- 
incide well with the positions of the northern and southern cores. 
If we interpret the velocity gradients as rotational motion, we see 
that the two cores are rotating in almost opposite directions. The 
angular velocity of the northern core is Q = 2.2 x 10"'^ rad 
s"' and the southern Q = 6.9 x 10"'^ rad s"' . ISteinacker et aP 
(120051) proposed that Oph D would have been formed by tur- 
bulent fragmentation. The directions and magnitudes of the gas 
velocity field derived from our observation support this view. 

A modified Bonnor-Ebert sphere approximates the southern 
core reasonably well. The density and temperature profiles are 
reproduced satisfactorily when the central number density of the 
model is He = 4 x 10* cm"^ and the model temperature drops 
from 10 K at the surface to 6. 1 K in the centre. Heavy ammonia 
depletion within 700 AU from the centre is required to succes- 
fully match the southern core model with the observation. 

The northern core is more challenging to model, and we did 
not find a single good solution. Two different parameter value 
combinations succeed in producing an acceptable match with the 
observation. The other one is a modified Bonnor-Ebert sphere 
with central density of = 2 x 10^ cm"^ and having ammonia 
abundance increase inwards, reaching X(NH3) = 2 x 10"** at 
the centre. An almost equally good match with the observation 
can be achieved with a more dense spheroidical core with tic - 
1x10* cm"^ and heavy ammonia depletion within the innermost 
1000 AU. Gas kinetic temperatures at the centres of the two core 
models are 8.8 K and 7.4 K, respectively. 

The locat ion of the cold spo t in the southern part of Oph D 
according to iHariu et aTl ( l2008h is almost in between our two 
cores. They observed the ground state (lio — > In) rotational 
transition of ortho-H2D^, which may peak at a different position 
than ammonia owing to chemistry. Nevertheless, the previous 
H2D'^ result is surprising as one might expect to find the coldest 
gas to concentrate in the core centres. 

The minimum line-of-sight average gas kinetic temperature 
derived from the interferometric NH3 observations is r^in - 
7.4 + 1.1 K. The difference between this result and the tem- 
perature reported by Hariu et al. (2008) (Tidn = 6.0 + 1.4 K) is 
not statistically significant, and can furthermore be understood 
by considering the following points. First of all, H2D^ probably 
traces denser, and thus colder gas than ammonia. Furthermore, 
our temperature represents the line-of-sight average, which in- 
evitably increases the temperature estimate as ammonia thrives 
also in warmer envelope of the core. The fact that we could not 
utilise the full resolution interferometric map, but were forced 



to use a map cleaned with a 10" restoring beam in our analysis, 
also influenced our results. 

The cores seem to be gravitationally bound. The observed 
angular velocities of both cores, especially the northern one, are 
significant. If there is enough time for the cores to evolve, the 
rotation may lead to a formation of a disk. The southern core 
simulation results suggest that the core may be eventually col- 
lapsing, as the dimensionless radius of the Bonnor-Ebert sphere 
is ^ = 6.75, slightly above the critical value. There is no ob- 
servational data to aid in the discussion regarding the role of 
magnetic field in the support of these cores. Some mechanisms, 
like a slight increase in the external pressure and, in the case of 
magnetic support, ambipolar diffusion is needed to trigger the 
collapse process. 
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